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Abstract Few- and many-fermion systems on the verge of stability, and consisting of strongly inter¬ 
acting particles, appear in many areas of physics. The theoretical modeling of such systems is a very 
difficult problem. In this work we present a theoretical framework that is based on the rigged Hilbert 
space formulation. The few-body problem is solved by exact diagonalization using a basis in which 
bound, resonant, and non-resonant scattering states are included on an equal footing. Current experi¬ 
ments with ultracold atoms offer a fascinating opportunity to study universal properties of few-body 
systems with a high degree of control over parameters such as the external trap geometry, the number 
of particles, and even the interaction strength. In particular, particles can be allowed to tunnel out of 
the trap by applying a magnetic-field gradient that effectively lowers the potential barrier. The result 
is a tunable open quantum system that allows detailed studies of the tunneling mechanism. In this 
Contribution we introduce our method and present results for the decay rate of two distinguishable 
fermions in a one-dimensional trap as a function of the interaction strength. We also study the nu¬ 
merical convergence. Many of these results have been previously published (R. Lundmark, C. Forssen, 
and J. Rotureau, arXiv: 1412.7175). However, in this Contribution we present several technical and 
numerical details of our approach for the first time. 
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1 Introduction 

The tunneling of particles, energetically confined by a potential barrier, is a fascinating quantum 
phenomenon which plays an important role in many physical systems. An exciting recent development 
in the context of multiparticle tunneling is the experimental realization of few-body Fermi systems 
with ultracold atoms [1^ ■ These setups are extremely versatile as they are associated with a high 

degree of experimental control over key parameters such as the number of particles and the shape of 
the confinirg potential. In addition, the interaction between particles can be tuned using Feshbach 
resonances 0, which in the case of trapped particles turns into a confinement-induced resonance [13 : 
Flflj l. The resulting interparticle interaction is of very short range compared to the size of the systems, and 
can be modeled with high accuracy by a zero-range potential. Such tunable open quantum systems 
provides a unique opportunity to investigate the mechanism of tunneling as a function of the trap 
geometry and the strength of the interparticle interaction. 
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Fig. 1 Panel (a): Trap potential, indicating the position of s.p. and two-body resonance states. Panel (b): 
Complex-momentum contour and Berggren basis states, highlighting the position of the s.p. resonance pole. 


In this Contribution we will consider a system of interacting, two-component fermions in a finite- 
depth potential trap. The trap is not deep enough to support a single-particle (s.p.) bound state, 
but does provide a quasi-bound state with a finite lifetime. We employ an effective ID potential 
corresponding to the stated potential for the experimental setup in Ref. (20| 
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where pVo, xr, and B' are tunable parameters and CB,a ~ 1- This potential is illustrated in Fig.[T](a). 

We recently introduced the rigged Hilbert space approach to the study of such open quantum 
systems Q . This method extends beyond the domain of Hermitian quantum mechanics and includes 
also time-asymmetric processes such as decays (see e.g. Ref. @ and references therein). In nuclear 
physics this formulation has been employed in the Gamow Shell Model BE BE Els to study 
threshold states and decay processes. Recently, it has also been used to model near-threshold, bound 
states of dipolar molecules Q. 


2 Method 

We will obtain solutions of the many-body Hamiltonian for interacting particles using an expansion 
of s.p. states in the so-called Berggren basis B- This complex-momentum basis includes S-matrix 
poles (bound and resonant states) as well as non-resonant scattering states. The use of this basis, 
constituting a rigged Hilbert space, is key to our approach as it allows to consistently include the 
continuum when finding eigensolutions of the open quantum system. The corresponding completeness 
relation is a generalization of the Newton completeness relation [T^ (defined only for real energy states) 
and reads 
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where |u„) correspond to poles of the S-matrix, and the integral of states along the contour L+, extend¬ 
ing below the resonance poles in the fourth quadrant of the complex-momentum plane, represents the 
contribution from the non-resonant scattering continuum [I| . This basis also features a non-conjugated 
inner product. 

Single-particle basis In order to generate the s.p. Berggren basis to be used in the many-body calcu¬ 
lation, we start with a complex-momentum basis 


{x\ k, n) 


sin(fca;) {II = 0) 
^ cos{kx) {n = 1 ) 


(3) 


that is equivalent to a combination of left- and right-travelling plane waves. In fact, this set corresponds 
to the Berggren basis for a s.p. Hamiltonian with H = 0. In practice, the momentum integral of 
Eq. (21) is performed using Gauss-Legendre quadrature and the momentum basis states are discretized 
accordingly. 

We will use the short-hand notation \ki) to denote an eigenstate with complex momentum ki and 


parity 77^. The matrix representation 


hij — {ki\h \ kj'^ 

(4) 

of the one-body Hamiltonian 


h{x) = —tT?/{2m)d^/dx^ -\- V{x) 

(5) 


is in general not Hermitian and not symmetric. The latter property can be recovered by redefining 
hij = where Wij are the Gauss-Legendre weights. The complex-momentum contour L'^ is 

selected such that the s.p. resonance energy, which is one of the eigenvalues, appears above the contour 
in the fourth quadrant. A specific example is shown in Fig. [H^b), in which the contour L'^ consists of 
four segments and is truncated at fc = fcmax- The s.p. resonance eigenstate |ures) can be described as 
a Gamow state Such a state behaves asymptotically as an outgoing wave with a complex-energy 
Cres = Cr — ijrl‘2- The imaginary part of the energy corresponds to the decay width 7 ^ and gives 
the half-life of the s.p. state, ti /2 = ln( 2 )fi/ 7 j., and the s.p. tunneling rate 71 = 7 r-//i. The full set of 
eigenvectors to the s.p. Hamiltonian includes the resonance state, as well as non-resonant scattering 
states with complex momenta very close to the original contour. Together, these eigenstates correspond 
to our s.p. basis states Ui = {|Mi)} [13 ■ 


Many-body problem The interaction between fermionic atoms in different hyperfine states is modeled 
by the zero-range potential Vi 2 (xi, X 2 ) = gb (xi — X 2 ), with g the tunable interaction strength. The 
fermions will be referred to according to their hyperfine spin state as “spin-up” (t) and “spin-down” 
( 4 ,), thus making an obvious connection with systems of spin 1/2 particles (e.g. electrons or nucleons). 

We will consider two particles in the trap, being the simplest example of a many-body system. The 
Hamiltonian is 


2 

H{xi,X2) = ^ 

i=l 




g5{xi 


X2)- 


( 6 ) 


The two-particle basis T 2 is then naturally constructed from the s.p. bases for the spin-up and -down 
fermions as T 2 = Gi(t) ® Gi(4,). Note that the spin-dependent CB^^-term in the trap potential (2]) will, 
in general, result in different s.p. bases for different spin states. In the following we will not directly 
compare to experimental results and will therefore restrict ourselves to spin-independent trap potentials 
with CB.a = 1- Let us first consider the situation of two non-interacting particles, i.e. 5 = 0. In this case, 
the ground state of the system, corresponds to the two distinguishable fermions both occupying 

the resonant (quasi-bound) state |Mres)- In this configuration, both particles are localized in the trap 
for a finite amount of time, before tunneling out through the potential barrier. 

In order to construct the many-body Hamiltonian matrix we first need to evaluate the interaction 
matrix elements between Berggren states. Note that the basis functions along the complex contour 
diverge for x —>■ 00 . As a consequence, the matrix elements of the two-body interaction are not finite 
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in the Berggren basis. We solve this issue by regularizing the two-body matrix elements between states 
in T 2 using an expansion in the harmonic oscillator (HO) basis Q 

'^max 

{ui,Uj\Vi2\uuUm) = ^ {Ui\na) {Uj\ni}) {ui\n^) {Urn\ns) {na,nj}\Vi 2 \n^,ns), (7) 

■) 

Tlry , 7Z(5 


where the sum runs over all HO states up to some truncation Umax- In the end, our Hamiltonian ([6]) 
matrix in this rigged Hilbert space will be non-Hermitian, but complex symmetric. The spectrum will 
include bound, resonant and scattering many-body states. 


Diagonalization We will now turn to the problem of finding the resonance state in the eigenspectrum 
of the many-body Hamiltonian matrix. While there are several algorithms available for finding extreme 
eigenvalues of Hermitian matrices, our problem is different. We are searching for many-body resonance 
solutions, l^res), that are characterized by outgoing boundary conditions and a complex energy iS^es = 
Er — irr/2, where Fr is the decay width due to the emission of particles out of the trap. In general, these 
physical states will correspond to specific complex eigenvalues in the interior part of the eigenvalue 
spectrum. Such eigenstates can be identified by the property that they will be independent of the 
particular choice of L+ as long as the Berggren completeness relation (O holds, i.e. fc^ax and the 
number of discretization points both need to be large enough. 

However, there exists a simpler method to distinguish these states from the continuum of many- 
body scattering solutions. The resonance state is usually the state with the largest overlap (in modulus) 
with 1^*-°^), referred to as the pole approximati on |lOl | . With the aim of targeting this state we employ 
the Davidson algorithm for diagonalization lll|. The Davidson method is very efficient at finding 
eigenvalues for diagonally dominant matrices. The basic idea of this algorithm is that a search space can 

be constructed by targeting a certain eigenpair. An approximation Afc j for the desired eigenpair 

E) is constructed in a subspace of dimension k that is much smaller than the full dimension. The 
search space is extended iteratively, as in many other methods, but the Davidson algorithm does not 

rely on a Krylov subspace. At each iteration, we select the Ritz pair that has the largest 

overlap with the pole approximation, |^(°^). The search space is then extended in the direction of the 
residual vector jr^) = E[\(j)k) — Ek\(j)k)- Convergence is achieved when the norm of the residual vector 
approaches zero. The main computational cost of the method is the matrix-vector multiplication that 
is required at each iteration. For the problems that we consider here, convergence is usually reached 
within 10-20 iterations as demonstrated in Fig. [21 


Many-body tunneling rate Concerning the tunneling rate we want to stress that there is a priori no 
simple relation between the decay width and the half life for a many-body system, contrary to the case 
of a s.p. Gamow state. Assuming exponential decay we would estimate the tunneling rate yr = /7-/h = 
—2Im(£’res)/Ii. However, having access to the resonant wave function, <l>T-es{xi,X 2 ) = ^res(x), we can 
alternatively compute the decay rate using an integral formalism Q. The rate of particle emissions 
can be obtained by integrating the outward flux of particles at large distance ccout from the center of 
the trap, and normalizing by the number of particles on the inside 


Tflux — 
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with A^(a;out) = Hi dxj\<Pres{^)\‘^■ 
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Fig. 2 Convergence of the resonance state for three different interaction strengths using the Davidson method. 
The model space dimension is 57600, corrsponding to A^pts = 240. Panel (a): Residual norm of the eigenstate. 
Panel (b): Relative convergence (distance from converged value) of the real part of the eigenvalue. Panel (c): 
Same as (b) but for the imaginary part of the eigenvalue. The converged eigenvalue is denoted kr,oo- 


3 Results 

In this Contribution we restrict ourselves to the simplest instance of the described tunable open 
quantum system, the case of two interacting fermions in different spin states in an open ID potential 
trap. However, we want to stress that the formalism can be applied to higher-dimensional traps and 
to systems with more particles. For comparison with experimental results we will use molecular units, 
in which energy is given in nKhn, time in /rs, and distances in fj,m. In these units we have h = 
7638.2 nKke ^s, the Bohr magneton = 6.7171 • lO^nKkeT”^ and f?jm = 80.645 nKke/rm^, 
where m is the mass of a ®Li atom. 

In Fig. [Ija) we show for illustrative purpose the trap potential with pVb = 2.123 • lO^nKhe, 
xb = 9.975 /rm, B' = 18.90 • 10“® T/rm“^, CB,a = 1, which closely resemble the parameters extracted 
from experimental data (see also discussion below). In order to handle the linear term B'x we truncate 
the potential at Xcut, sufficiently far away from the relevant trap region. In practice, this is achieved 
by applying a positive energy shift Fishift so that ^(xcut) + F'shift = 0. The energy shift is subtracted at 
the end, and we have verified that the fluctuations in the s.p. energy (tunneling rate) with the choice 
of Fishift was less than 0.04 % (2%). 

The s.p. Schrbdinger equation is solved using the method described above. The discrete set of 
complex-momentum states {|fci)} that span the contour F+ is shown as blue dots in Fig. [Ijb). The 
energy shift that was used is Bghift = bOOuKke- The resulting set of eigenstates (green circles) lies 
very close to the contour with the exception of one isolated state. The former states correspond to 
non-resonant scattering solutions, while the latter is a resonance. Together, these eigenstates form the 
complete set of s.p. basis states, {|ui)}, that will be used in the many-body calculation. 
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Fig. 3 Two-fermion resonance state as a function of the interaction strength g for = 1. Panel (a): 
Interaction energy (0 compared with the corresponding energy obtained using the pole approximation. Panel 
(b): Tunneling rates obtained from the imaginary part of the resonance energy (from the full calculation and 
the pole approximation, respectively) compared with the rate obtained from the flux calculation (jS)). 


The number of points on the contour is increased until convergence of the s.p. resonance energy 
is achieved. Note that the resonance pole will always remain fixed while the set of scattering states 
will depend on the choice of the contour L+. For illustrative purposes, the contour shown in Fig. [1] 
consists of only TVpts = 100 basis states while full calculations were performed with iVpts = 240-320. 
For this set of potential parameters we hnd Cres = (301.415 — 0.08548z) uKke, which translates into a 
tunneling rate 71 = 22.38 s“^. 

We now consider the solution of the interacting two-fermion system, projected on the full Berggren 
basis. We define the interaction energy as 

Eint = Re(£'res) ~ 26^, (9) 

where Re(ifres) = E^ is the real part of the resonance energy. Results for the two-particle resonance 
state as a function of the interaction strength g are shown in Fig. [3] For g = 0, the two fermions tunnel 
out independently and the tunneling rate is equal to 7 = 271 = 44.76 s~^. However, as the interaction 
becomes more attractive, the real part of the resonance energy decreases, and the effective barrier seen 
by the two particles increases. As a consequence the tunneling rate decreases as seen in Fig. [3Kb). 

Along with the full calculations, we show in Fig.|3]also results obtained in the pole approximation, 
which corresponds to the single configuration where the two distinguishable fermions occupy the s.p. 
resonant state. This comparison clearly demonstrates the importance of continuum correlations. The 
resonance energy and width are both decreased due to configuration mixing between the s.p. resonance 
pole and non-resonant scattering states. In particular, the energy width, which translates into a decay 
rate, is very sensitive to these correlations. These results highlight the importance of properly taking 
the openness of the system into account. 

The agreement between the tunneling rate computed from the decay width of the resonance and 
from the flux formula ([3]) demonstrates the quality of our numerical approach. It also shows that the 
tunneling is well approximated by an exponential decay law for this system. 


4 Numerical convergence 

The stability of the results for the resonance energy, with respect to different model-space param¬ 
eters, may be investigated in order to assess the numerical precision of the method. We have per¬ 
formed a series of such convergence studies for systems with different interaction strengths {g = 
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Fig. 4 Convergence stndy for the two-body problem with different interaction strengths. The upper row shows 
relative changes of the resonance pole position with different contours and different number of dicretization 
points, for strong attraction (panel a) and intermediate attraction (panel b). The two lower rows show the 
relative change of the real part (panel c) and imaginary part (panel d) of the resonance energy with increasing 
HO trnncation in the calculation of two-body matrix elements 0. 


-blOO, —20, —100 nK kB/im). Panels (a) and (b) of Fig. H] show the relative change of the resonance 
energy position for different contours L'^ and different number of discretization points. Panels (c) and 
(d) of Fig. 0] demonstrate the convergence with increasing HO truncation in the computation of inter¬ 
action matrix elements ©• Based on these results we quantify the numerical precision for a specific 
interaction strength by adding (in quadrature) the amplitudes of variations when the model-space pa¬ 
rameters were altered one by one. We found that the precision of the real part of the interaction energy 
was on the order of < 2 % for the entire range of interaction strengths. However, the precision of the 
computed imaginary energy was found to have a lower bound since variations of the computed decay 
rate was never smaller than 0.5 s“^. This becomes obvious when the interaction is strongly attractive 
and the absolute value of the decay rate is (~1 see Fig. SKa). 

The larger relative variations in the tunneling rates compared to the interaction energies, as demon¬ 
strated in this section, are most likely due to the fact that the imaginary part of the energy is several 
orders of magnitude smaller than the real part. This means that in order to get precise results for the 
tunneling rates, an even more precise result for the modulus of the energy is needed. The estimated 
uncertainties from these numerical studies are shown as shaded bands in both panels of Fig. [31 but is 
only visible in the tunneling rate for the most attractive interactions {g < — GOnKke /im). 


5 Concluding remarks 

The tunneling of few fermions from low-dimensional traps were measured by Ziirn et al [2^ . The anal¬ 
ysis of data from this experiment is quite complicated and involves the use of the WKB approximation 
to extract the trap potential parameters. More precisely, pVo and B' in Eq. 0 were adjusted such 
that the s.p. tunneling rates obtained in the WKB approximation matched the experimental results. 
We studied this experiment in a recent publication Q using the method outlined in this Contribution. 
Using the set of parameters given in Ref. as input to our exact diagonalization approach we found 
a good agreement for the s.p. energies (with a difference of at most a f ew p ercent), while s.p. tunneling 
rates were almost two times larger than the ones published in Ref. |20| . The main conclusion from 
these hndings is that that the WKB method should not be expected to produce reliable estimates for 










the tunneling rate and that the analysis of experimental results for open quantum systems is highly 
sensitive to the determination of trap parameters. 


Acknowledgements The first author wishes to acknowledge the organizers and participiants of the 7th 
International Workshop on the “Dynamics of Critically Stable Quantum Few-Body Systems” in Santos, Brazil. 
The cross-disciplinary theme of the meeting provided a very stimulating environment and triggered many 
interesting discussions. The research leading to these results has received funding from the European Research 
Council under the European Community’s Seventh Framework Programme (FP7/2007-2013) / ERC grant 
agreement no. 240603, and the Swedish Foundation for International Cooperation in Research and Higher 
Education (STINT, IG2012-5158). The computations were performed on resources provided by the Swedish 
National Infrastructure for Computing (SNIC) at High-Performance Computing Center North (HPC2N) and 
at Chalmers Centre for Computational Science and Engineering (C3SE). We thank the European Centre for 
Theoretical Studies in Nuclear physics and Related Areas (ECT*) in Trento, and the Institute for Nuclear 
Theory at the University of Washington, for their hospitality and partial support during the completion of this 
work. We are much indebted to D. Blume, M. Zhukov, N. Zinner, and G. Ziirn for stimulating discussions. 


References 

1. Berggren T (1968) On the use of resonant states in eigenfunction expansions of scattering and reaction 
amplitudes. Nucl Phys A109:265-287 

2. Chin C, Grimm R, Julienne P, Tiesinga E (2010) Feshbach resonances in ultracold gases. Rev Mod Phys 
82(2):1225-1286 

3. Davidson ER (1975) Iterative Calculation of a Few of Lowest Eigenvalues and Corresponding Eigenvectors 
of Large Real-Symmetric Matrices. J Comput Phys 17(l):87-94 

4. Fossez K, Michel N, Nazarewicz W, Ploszajczak M (2013) Bound states of dipolar molecules studied with 
the Berggren expansion method. Phys Rev A 87(4):042,515 

5. Gamow G (1928) Zur Quantentheorie des Atomkernes. Zeitschrift fiir Physik 51(3):204-212 

6. Grigorenko LV, Zhukov MV (2007) Two-proton radioactivity and three-body decay. HI. Integral formulas 
for decay widths in a simplified semianalytical approach. Phys Rev C 76(1):014,008 

7. Hagen G, Hjorth-Jensen M, Michel N (2006) Gamow shell model and realistic nucleon-nucleon interactions. 
Phys Rev C 73(6):064,307 

8. Lundmark R, Forssen G, Rotureau J (2014) Tunneling Theory for Tunable Open Quantum Systems of 
Ultracold Atoms in One-Dimensional Traps. ArXiv e-prints 1412.7175 

9. de la Madrid R (2005) The role of the rigged Hilbert space in quantum mechanics. Eur J Phys 26(2):287-312 

10. Michel N, Nazarewicz W, Ploszajczak M, Bennaceur K (2002) Gamow Shell Model Description of Neutron- 
Rich Nuclei. Phys Rev Lett 89(4):042,502 

11. Michel N, Nazarewicz W, Ploszajczak M, Okolowicz J (2003) Gamow shell model description of weakly 
bound nuclei and unbound nuclear states. Phys Rev C 67(5) 

12. Michel N, Nazarewicz W, Ploszajczak M, Vertse T (2009) Shell model in the complex energy plane. J Phys 
G 36(1):013,101 

13. Newton RG (1960) Analytic Properties of Radial Wave Functions. J Math Phys 1:319-347 

14. Olshanii M (1998) Atomic Scattering in the Presence of an External Gonfinement and a Gas of Impenetrable 
Bosons. Phys Rev Lett 81(5):938-941 

15. Papadimitriou G, Kruppa AT, Michel N, Nazarewicz W, Ploszajczak M, Rotureau J (2011) Gharge radii 
and neutron correlations in helium halo nuclei. Phys Rev C 84(5):051,304 

16. Rotureau J, van Kolck U (2013) Effective Field Theory and the Gamow Shell Model. The ®He Halo Nucleus. 
Few-Body Syst 54(5):725-735 

17. Rotureau J, Michel N, Nazarewicz W, Ploszajczak M, Dukelsky J (2006) Density Matrix Renormalization 
Group Approach for Many-Body Open Quantum Systems. Phys Rev Lett 97(11):110,603 

18. Serwane F, Ziirn G, Lompe T, Ottenstein TB, Wenz AN, Jochim S (2011) Deterministic Preparation of a 
Tunable Few-Fermion System. Science 332(6027):336-338 

19. Ziirn G, Lompe T, Wenz AN, Jochim S, Julienne PS, Hutson JM (2013) Precise Characterization of Li-6 
Feshbach Resonances Using Trap-Sideband-Resolved RF Spectroscopy of Weakly Bound Molecules. Phys 
Rev Lett 110(13):135,301 

20. Ziirn G, Wenz AN, Murmann S, Bergschneider A, Lompe T, Jochim S (2013) Pairing in Few-Fermion 
Systems with Attractive Interactions. Phys Rev Lett 111(17):175,302 



